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We present a set of three-dimensional (3D) direct numerical simulations of incompressible decaying 
magnetohydrodynamic turbulence in which we investigate the influence of an external uniform 
magnetic field Bo- A parametric study in terms of Bo intensity is made where, in particular, 
we distinguish the shear- from the pseudo- Alfven waves dynamics. The initial kinetic and magnetic 
energies are equal with a negligible cross-correlation. Both the temporal and spectral effects of 
Bo are discussed. A sub-critical balance is found between the Alfven and nonlinear times with 
both a global and a spectral definition. The nonlinear dynamics of strongly magnetized flows is 
characterized by a different /c^-spectrum (where Bo defines the parallel direction) if it is plotted at 
a fixed k\\ (2D spectrum) or if it is integrated (averaged) over all k\\ (ID spectrum). In the former 
case a much wider inertial range is found with a steep power law, closer to the wave turbulence 
prediction than the Kolmogorov one like in the latter case. It is believed that the averaging effect 
may be a source of difficulty to detect the transition towards wave turbulence in natural plasmas. 
For the first time, the formation of filaments is reported within current and vorticity sheets in 
strongly magnetized flows which modifies our classical picture of dissipative sheets in conductive 
flows. 

PACS numbers: 47.27.Jv, 47.65.-d, 52.30.Cv, 95.30.Qd 



I. INTRODUCTION 

The magnetohydrodynamics (MHD) approximation 
has proved to be quite successful in the study of a variety 
of astrophysical plasmas, electrically conducting gas or 
fluids, like those found in the solar corona, the interplan- 
etary medium or in the interstellar clouds. These me- 
dia are characterized by extremely large Reynolds num- 
bers (up to 10 13 ) d with a range of available scales from 
10 18 m to few meters. The isotropy assumption, usually 
used in hydrodynamic turbulence, is particularly difficult 
to justify when dealing with astrophysical flows since a 
large-scale magnetic field is almost always present like 
in the inner interplanetary medium where the magnetic 
field lines form an Archimedean spiral near the equato- 
rial plane (see e.g. [E @]). Thus, MHD turbulence is 
much more complex than Navier-Stokes turbulence with, 
in particular, a nonlinear transfer between structures of 
various sizes due to both nonlinear couplings and Alfven 
wave propagation along the background magnetic field. 

In the mid sixties, Iroshnikov Q and Kraichnan [5] 
(hereafter IK) proposed a first description of incompress- 
ible MHD turbulence. In this approach a la Kolmogorov, 
the large-scale magnetic field is supposed to act on small- 
scales as a uniform magnetic field, leading to counter- 
propagating Alfven waves whose interactions with turbu- 
lent motions produce a slowdown of the nonlinear energy 
cascade. The typical transfer time through the scales is 
then estimated as t^ l /ta (instead of tnl for Navier- 
Stokes turbulence), where tnl ~ tjui is the nonlinear 
eddy turnover time at characteristic length scale i and U£ 
is the associated velocity. The Alfven time is ta ~ £/Bq 
where Bo represents the large-scale magnetic field nor- 
malized to a velocity (Bo —> BoyT^oPo? with /io the 



magnetic permeability of free space and po the uniform 
plasma density). Note that we will use this renormal- 
ization in the rest of the paper. Hence, the IK energy 
spectrum in k~ 3 / 2 unlike the /c -5 / 3 Kolmogorov one for 
neutral flows. 

The weakness of the IK's phenomenology is the appar- 
ent contradiction between the presence of Alfven waves 
and the absence of an external uniform magnetic field. 
The external field is supposed to be played by the large- 
scale magnetic field but its main effect, i.e. anisotropy, 
is not included in the description. The role of a uni- 
form magnetic field has been widely discussed in the lit- 
erature and, in particular, during the last two decades 
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At strong B intensity, one of the most clearly estab- 
lished results is the bi-dimensionalization of MHD tur- 
bulent flows with a strong reduction of nonlinear trans- 
fers along Bq. In the early eighties, it was shown that 
a strong Bo leads to anisotropic turbulence with an en- 
ergy concentration near the plane k • Bo = [6], a re- 
sult confirmed later on by direct numerical simulations 
in two and three space dimensions 0, Q. A linear de- 
pendence between anisotropy and Bq intensity was also 
suggested [13) . From an observational point of view, we 
have also several evidences that astrophysical (and lab- 
oratory) plasmas are mostly in anisotropic states like in 
the solar wind (see e.g. 0, (24|) or in the interstellar 
medium (see e.g. p5|). 

The effects of a strong uniform magnetic field may be 
handled through an analysis of resonant triadic interac- 
tions [7] between the wavevectors (k, p, q) which satisfy 
the relation k = p + q, whereas the associated wave fre- 
quencies satisfy, for example, 00 (k) = cj(p) — 00 (q). The 
Alfven frequency is o;(k) = k B = k\\Bo, where || de- 
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fines the direction along B (_L will be the perpendicular 
direction to B ). The solution of these three- wave res- 
onant conditions directly gives, q\\ = 0, which implies a 
spectral transfer only in the perpendicular direction. For 
a strength of Bo well above the r.m.s. level of the kinetic 
and magnetic fluctuations, the nonlinear interactions of 
Alfven waves may dominate the dynamics of the MHD 
flow leading to the regime of (weak) wave turbulence 
where the energy transfer, stemming from three- wave 
resonant interactions, can only increase the perpendic- 
ular component of the wave vectors, while the nonlinear 
transfers is completely inhibited along B [3 EH . 

Another important issue discussed in the literature 
is the relationship between perpendicular and parallel 
scales in anisotropic MHD turbulence (see [1, [Io|, IZoj). 
In order to take into account the anisotropy, Goldreich 
and Shridar [l0[ proposed a heuristic model based on a 
critical balance between linear wave periods and nonlin- 
ear turnover time scales, respectively ta ~ £\\/Bq and 
tnl ~ £±/v>£ (where £\\ and £± are the typical length 
scales parallel and perpendicular to Bo), with ta = tnl 
at all inert ial scales. Following the Kolmogorov argu- 
ments, one ends up with a E(k_\_,k\\) ~ kj 6 ^ 3 energy 
spectrum (where k = (kj_,fci|) and k± = | kj_ | ) with the 
anisotropic scaling law 

h~k 2 '\ (i) 

A generalization of this result has been proposed re- 
cently [26[ in an attempt to model MHD flows in both 
the weak and strong turbulent regimes, as well as in 
the transition between them. In this heuristic model, 
the time-scale ratio x — t a/tnl is supposed to be con- 
stant at all scales but not necessarily equal to unity. The 
relaxation of this constraint enables to still recover the 
anisotropic scaling law (pQ) and to find a universal predic- 
tion for the total energy spectrum E(k±, fey) ~ kj_ a k^ , 
with 3a + 2(3 = 7. Accor ding to direct numerical sim- 
ulations (see, e.g. [2?], HH, l29j]), one of the most funda- 
mental results seems to be the anisotropic scaling law 
between parallel and perpendicular scales (pp) and an ap- 
proximately constant ratio x-> generally smaller than one, 
between the Alfven and the nonlinear times. This sub- 
critical value of x implies therefore a dynamics mainly 
driven by Alfven waves interactions. 

In the weak turbulence limit, the time-scale separa- 
tion, x < 1, leads to the destruction of some nonlinear 
terms, including the fourth-order cumulants, and only 
the resonance terms survive [3, EH, OiO, [3l| which allows 
to obtain a natural asymptotic closure for the wave ki- 
netic equations. In absence of helicities and for k± ^> fey , 
the dynamics is then entirely governed by shear- Alfven 
waves, the pseudo- Alfven waves being passively advected 
by the previous one. In the case of an axisymmetric tur- 
bulence, and in the absence of cross-correlation between 
velocity and magnetic field fluctuations, the exact power 
law solution is E(k±, fen) ~ kj 2 f(k\\), where / is an arbi- 
trary function taking into account the transfer inhibition 



along B . The regime of wave turbulence is quiet diffi- 
cult to reproduce by direct numerical simulations since it 
requires a strong external magnetic field as well as a high 
spatial resolution. According to a recent theoretical anal- 
ysis, it seems to be currently not possible to reach fully 
this regime [32| and only the transition towards such a 
regime is likely to be obtained [22|, [33[ . 

In order to better understand the development of 
anisotropy in natural magnetized plasmas, we perform 
a set of tri-dimensionnal numerical simulations of incom- 
pressible MHD. In this work, the regime of freely decay- 
ing flows is chosen in an attempt to model the nonlin- 
ear evolution of outward and inward propagating Alfven 
waves. We mainly focus our analysis on the develop- 
ment of anisotropy in flows at moderate Reynolds num- 
bers which freely evolve under the influence of a uniform 
magnetic field whose strength will be taken as a param- 
eter. The details of the numerical setup and simulations 
are given in the next Section. In Section III, we inves- 
tigate the temporal characteristics of the different flows, 
as well as different global quantities to measure the spec- 
tral anisotropy. Section IV is devoted to the evolution of 
the energy spectra together with their fluxes. The flow 
spatial properties are examined in Section V. Section VI 
discusses the second set of simulations. A summary and 
a conclusion are finally given in Section VII. 

II. NUMERICAL SETUP 

A. Incompressible MHD equations 

The MHD equations that describe the large-scale and 
low frequency dynamics of magnetized plasmas are, in 
the incompressible case and in the presence of a uniform 
magnetic field B , 

d t v - B d {l b + v ■ Vv = -VP* + b- Vb + z/Av , (2) 
d t b - B d\\v + v- Vb = b- Vv + rjAb , (3) 
V • v = , (4) 

Vb = 0, (5) 

where v is the plasma flow velocity, b the magnetic field 
(normalized to a velocity), P* the total (magnetic plus 
kinetic) pressure, v the viscosity and r\ the magnetic dif- 
fusivity. It is convenient to introduce the Elsasser fields 
z ± = u±b for the fluctuations; in this case and assuming 
a unit magnetic Prandtl number (i.e. ^ = r/), we get 

d t z ± + z T • Vz ± =f P S||Z ± = - VP* + ^V 2 z ± , (6) 
V-z ± = 0. (7) 
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Note that the second term in the left hand side (LHS) 
of Eq. (j6]) represents the nonlinear interactions between 
the z ± fields, while the third term represents the linear 
Alfvenic wave propagation along the Bq field which will 
be assimilated to the z-direction in our numerical box. 
In the present analysis, a unit magnetic Prandtl num- 
ber is taken in order to extend at maximum the inertial 
range for both the kinetic and magnetic energies. We 
believe that such analysis is the first step to understand 
turbulence in anisotropic media. The extension to other 
magnetic Prandtl numbers is the second step: this situa- 
tion, more realistic for turbulence like in the interstellar 
medium, supposes to keep a high level of turbulence for 
both the kinetic and magnetic energies which necessitates 
a higher spatial resolution. 



1. Runs la to IVa 

The initial kinetic and magnetic fluctuations are char- 
acterized by spectra at large-scales, i.e. for k = [1,8], 
proportional to k 2 exp(— /c 2 /4); for k > 8, the spectra 
are exactly equal to zero. This condition means that for 
wavenumbers k up to 2, we have mainly a flat modal 
spectrum which prevents initially any favored wave vec- 
tors. No forcing is present during the simulations and the 
flows may evolve freely for time t > 0. The associated 
kinetic, 



£™ = -<u 2 (x)>, 



(11) 



and magnetic, 



B. Poloidal/toroidal decomposition 

In the presence of a large-scale magnetic field Bo, 
Alfven waves develop and propagate at Alfven speed Bq 
along the Bq direction. These waves may be decom- 
posed into shear- and pseudo- Alfven waves denoted, re- 
spectively, zf and • The divergence- free condition im- 
plies that only two types of scalar field and ± ) are 
needed to describe the incompressible MHD dynamics 
which are, respectively, the toroidal and poloidal fields. 
For the Fourier transforms of the involved fields, we have: 



: (k) 



2±(k) 



with 



zf(k) = ik x e||'^ ± (k) , 



: (k) = 



k x (k x en) . 



: (k). 



(8) 

(9) 
(10) 



where in our simulations k± = A / k 2 + k 2 and k 



kx + ky + k 2 . Here, k z = k\\ and ey denotes the unit 

vector parallel to the B direction. Hence, the shear- 
Alfven waves correspond to a vector field perpendicular 
to the external magnetic field B whereas the pseudo- 
Alfven waves is a vector field which may have a compo- 
nent along Bq; but both vector fields depend on the three 
coordinates of k. 



C. Initial conditions 

We numerically integrate the three-dimensional incom- 
pressible MHD equations ©-(O, in a 27r-periodic box, 
using a pseudo-spectral code (including de-aliasing), and 
with spatial resolution from 256 3 to 51 2 2 x 64 grid points 
according to the initial conditions (see Table [Ij). The time 
marching uses an Adams-Bashforth / Cranck-Nicholson 
scheme, i.e. a second-order finite-difference scheme in 
time (see e.g. [HI). 



£ b = -<b 2 (x)>, 



(12) 



energies are chosen initially equal, namely E v (t = 0) = 
E b (t = 0) = 0.5. (Note that < • > means space averag- 
ing.) 

The correlation between the velocity and magnetic 
field fluctuations, which is measured by the cross- 
correlation 



. 2 < u(x) • b(x) > 
= <u 2 (x)+b 2 (x) > ' 



(13) 



is initially less than 1%. 

The initial (large-scale) kinetic and magnetic Reynolds 
numbers are about 800 for the flows with v = 4 x 10 -3 



(see Table % 
gral scale is 



with u r 



1; the isotropic inte- 
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f(Ev(k)/k)dk 
fE-(k)dk 



(14) 



A parametric study is performed according to the in- 
tensity of Bo. Four different values are used, namely 
Bq = 0, 1, 5 and 15. All these simulations are run up to a 
maximum computational time Im — 15, and correspond 
to runs la to IVa described in Table H 



2. Runs Va to Vila 

Taking advantage of the strong reduction of the non- 
linear transfers along Bq in highly magnetized flows, a 
second set of direct numerical simulations is performed 
with a spatial resolution of 512 2 grid points in the per- 
pendicular plane to B and with only 64 grid points in 
the parallel direction (runs Va and Via in Table [Ij). For 
such runs, the initial conditions are the same as before 
with, however, a uniform magnetic field Bo = 15 and 30, 
and a smaller viscosity. 

Such simulations were analyzed in the past to explore 
the self-consistency of the reduced MHD model [35[ with 
the conclusion that small values of viscosities, adjusted 



TABLE I: Computational parameters for runs la to Vila with isotropic initial conditions, and for runs lb and lib with specific 
initial conditions (see text). Note that simulations Vila and lib use a hyperviscosity and a hyperdiffusivity (dissipation terms 
in V 4 ). Spatial resolution, viscosity v (= 77) and applied magnetic field intensity Bo are given, followed by initial integral length 
scales: isotropic L = 2tt f (E v (k) /k)dk/ f E v (k)dk, perpendicular L± — 2tt f (E v (k±) /k±)dk±/ f E v (k±)dk±, and parallel 
L\\ = 2tt f (E v (k\\) /k\\)dk\\/ f E v (k\\)dk\\ scales. Initial r.m.s. velocity u rms =< v 2 > 1//2 (= b rms =< b 2 > 1 / 2 ) fluctuation is 
given together with the initial kinetic Reynolds number 1Z V = u r m S L/v. Finally, we find typical times: isotropic eddy turnover 
time t nl = LI u rrns (based on the isotropic length- scale L), eddy turnover time tnl = L±/u r ms (based on Alfven time 

based on r.m.s. magnetic fluctuations r\ — L/b r ms, Alfven wave period ta — L\\/Bq and the final time £m of the numerical 
simulation. 
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according to the transverse dynamics, are not incompat- 
ible with the smaller spatial resolution in the parallel 
direction since the transfer toward small-scales is also re- 
duced along the uniform magnetic field. We checked that 
the viscosity, v = 10 -3 , is indeed well adjusted. Note 
that such a small aspect ratio may reduce the number of 
resonant wave interactions which in turns may affect the 
dynamics (see e.g. (36[). However, in Alfven wave tur- 
bulence, the resonant manifolds foliate wave vector space 
[lij which, in principle, prevent such a problem. 

In the same manner, another computation (run Vila) 
is made using an hyperviscous scheme, where the lapla- 
cian operator of the dissipative terms is replaced by a 
bi-laplacian, in order to enlarge the inertial range of the 
energy spectra. 



3. Runs lb and lib 

Finally, to evaluate the influence of the initial condi- 
tions, a third set of runs is performed with a uniform 
magnetic field fixed to Bq = 15, and with either a vis- 
cous (lb) or an hyperviscous dissipation (lib). In both 
simulations, we use a 512 2 x 64 grid points. The specific 
initial conditions of these runs correspond to a modal 
energy spectrum E ± (k±, fey ) = C(k\\ )/c 3 L , for k± and fey 
G [0,4], the value of C(fcy) increasing with fey to reach 
a maximum at fey = 4. Note that this initial spectrum 
allows a transient period of cascade toward smaller scales 
during which energy is mainly conserved. Initially, the 
ratio between kinetic and magnetic energies is still fixed 
to 1, whereas the cross-correlation coefficient is zero. A 
first set of results was given in [22I ]. 

For all the runs described in this Section, the com- 
putational parameters (initial Reynolds numbers, char- 
acteristic length scales and times...) are summarized in 
Table H 



III. TEMPORAL ANALYSIS 
A. Energetic properties 

1. Elsasser z ± cartesian fields 

In this section, we study the temporal behavior of 
several global quantities to characterize the MHD flow 
dynamics and the influence of the Bq strength on it. 
In all the following figures, time evolutions are shown 
from initial isotropic conditions up to time tu (the max- 
imum computational time reached), for simulations Ia 
to IVa at moderate resolution (256 3 mesh points) and 
Bo = 0,1,5 and 15, together with highly magnetized 
flows Va and Via at Bo = 15 and 30, using 512 2 x 64 
spatial resolution (see Table Ij). 

We first consider the evolution of the Elsasser energies, 

i? ± (t) = i<z ±2 (x)>(t), (15) 

displayed in Figure [TJ Note that, for periodic boundary 
conditions, these energies are two independent invariants 
of the inviscid MHD equations (|6]), with or without the 
presence of a uniform magnetic field. Energies E + (t) and 
E~(t) present a similar behavior for a given Bq. For runs 
Ia to Ilia, where B intensity is increased, we clearly 
see a slowdown of the energy decay. On one hand, this 
slowing down reflects the energy transfer inhibition along 
the Bq direction, and thus, the flow inability to create, 
in the parallel direction, smaller and smaller scales up 
to the dissipative ones. Hence, the energy dissipation 
mainly takes place in transverse planes which are led to 
play a more efficient role as the flow magnetization is in- 
creased. On the other hand, energy transfers themselves 
could also be weakened (in the transverse planes) since 
the MHD cascade of energy to smaller scales is produced 
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FIG. 1: Temporal evolution of energies E + (top) and E 
(bottom) for B = 0, 1, 5, 15 (runs la to IVa; 256 3 ) and B = 
15,30 (runs Va and Via; 512 2 x 64). The hyperviscous run 
Vila with Bq = 15 (512 2 x 64) is also given up to t = 14. 



by successive interactions of oppositely directed waves. 
Indeed, for higher Bq intensities, the waves become faster 
and thus the time duration of individual collision of z ± 
waves decreases. Therefore it takes much more collisions 
between (fast) Alfven wave packets (as measured by the 
ratio between the nonlinear turnover time on the linear 
wave period; tnl/ta) to have an efficient energy cascade 
process. One could also note that, for a given flow, a 
saturation effect occurs according to Bo intensities. In- 
deed, the E ± (t) evolutions are quite similar for flows at 
v = 4.10 -3 with Bo = 5 and 15 (runs Ila and IVa, re- 
spectively), as well as for flows at v = 10 -3 with Bo — 15 
and 30 (runs Va and Via, respectively). The hypervis- 
cous run Vila is also shown but only up to t = 14. We 
see that the initial plateau is wider and almost flat be- 



FIG. 2: Temporal evolution of the global dissipation vQ^ 
(top) and (bottom); same viscous runs as in Figure [1] 



cause of the larger inertial range and the higher Reynolds 
number. Then, we see a decay of energy which is slower 
than for the other viscous runs. 

The Bo saturation effect is also visible on the time 
evolution of the global dissipation of the flow, 



z/fi ± (t) = v < [V x z±] 2 (x) > (t) , 



(16) 



displayed in Figure [2j The early time dynamics, near the 
first inflection point, is almost inviscid; it corresponds to 
the small-scale generation (e.g. at times t < 1 in the 
Bo = simulation). As the Bq intensity is increased, 
this small-scale development is slightly retarded which 
means that the duration of the essentially inviscid phase 
increases. Moreover, the maximum of the dissipation is 
substantially reduced, and occurs at later times, namely 
t ~ 2 for flows with Bo = and 1, t ~ 3 with Bo = 5, 
15 and t ~ 4 for the less viscous flows with Bq = 15, 
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FIG. 3: Temporal evolution of the cross-correlation coefficient 
p\ same runs as in Figure [T] 



30. Altogether, in physical space, this corresponds to the 
creation of more elongated structures along B as the 
flow is more magnetized, with a smaller dissipation on 
the whole, and, in spectral space, to higher inhibition of 
parallel energy transfers, as already explained. One can 
also note that the dissipation peak is smoothed in the less 
viscous flows (y = 10 -3 ), meaning an almost constant 
dissipation between t ~ 3 and t ~ 5 in runs Va and Via 
with a more extended range of small scales. Finally, note 
a different evolution between case IVa (with v = 4.10 -3 ) 
and Va (with v = 10 -3 ) whereas the uniform field Bq 
is the same. A factor 4 of difference is visible initially 
which may be attributed mainly to a decrease of factor 4 
of the viscosity. In this case, the time delay to reach the 
maximum may be explained by a wider inertial range 
in k± and therefore a longer time needed to reach the 
dissipative scales (an effect also seen in Figure Q] with a 
wider initial plateau where energy is roughly conserved). 

Figure [3] shows the cross-correlation coefficient ([T3]) be- 
tween velocity and magnetic fields which also reads, in 
terms of the Elsasser energies, as 



Pit) 



E + (t)-E-(t) 
E+(t)+E-(t) ' 



(17) 



It measures the relative amount of the two z ± species. 
Indeed, p(t) — > ±1 means that E T = 0, and hence only 
one type of waves is excited, whereas when p(t) — > 0, they 
are as many z + as z~ counterpropagating waves, with the 
same amount of energy. Initially, pit = 0) ~ (i.e. less 
than 1%), and stays so during the flow inviscid phases. 
Close to the times at which the maximum of dissipation 
occurs in the different flows, p(t) deviates from zero with 
a lesser departure as the flow is more magnetized, from 
Bq = 1 up to Bq = 30, because the field lines are rigidi- 



fied by the ambiant magnetic field, and the dissipation is 
delayed. In the case of Bq =0, the temporal evolution of 
the cross-correlation coefficient is globally different, due 
to the absence of a guiding magnetic field and different 
dissipative processes. Note however that all flows evolve 
toward an excess of E~ energy. 
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4: Time evolution of the Alfven ratio ta for runs la to 



Apart from this, the prevalence of the Alfven wave fluc- 
tuations can be measured by the so-called Alfven ratio 



r A (t) = 



E v (t) _ < v 2 > (t) 

~mtj ~ < b 2 > jt) 



(18) 



between kinetic and magnetic energies. For example in 
the wave turbulence regime we have an equipartition (at 
the level of the kinematics [l4| ) between kinetic and mag- 
netic energies. Its departure from unity suggests the pres- 
ence of non Alfvenic fluctuations. Indeed, the energy of 
an individual Alfven wave is equipartitioned between its 
kinetic and magnetic components, averaged over a wave 
period, with thus a ratio ta = 1. In presence of an ex- 
ternal magnetic field, exchanges between magnetic and 
velocity fluctuations, due to Alfven waves, produce os- 
cillations as shown on r a (t) in Figure [H The period of 
these oscillations is given by the Alfven time ta ~ 5, 
1 and 0.4 (see Table [I]) which are found by a simple 
analysis based on the values Bq = 1, 5 and 15, respec- 
tively, and the values of the characteristic parallel length 
scale ZJ ~ 5.57 for runs Ila to IVa. Although, ini- 
tially the magnetic and kinetic energies are chosen equal, 
E v (t = 0) = E b (t = 0) = 0.5, the magnetic energy stabi- 
lizes around twice the kinetic energy, after t ~ 5, for the 
non-magnetized flow (Bq = 0). While for the magnetized 
flows, whatever the Bq intensity is, the magnetic energy 
saturates to about 1.25 lower than the kinetic energy 
level after time t ~ 2. This result may be compared with 
solar wind data where the same tendency is found with 
a domination of the magnetic energy. (This comparison 
is however not direct since outward propagating Alfven 
waves are initially dominant.) This Alfven ratio seems 
to find a limit of about 1 /2 at several astronomical units 
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which might be explained by the decreasing importance 
of the large-scale magnetic field at larger heliocentric dis- 
tances (see e.g. (37^). 

Figure [5] displays the pdf of the cross-correlation for 
different times (same runs as in Figured]). As expected, 
we start initially with a distribution clearly centered 
around zero. As the time increases, we see a distribu- 
tion shifted towards negative values to finally be centered 
around —0.4 for the non-magnetic case. The case Bq = 1 
is even more shifted with a maximum of the distribution 
around —0.6. The strongly magnetic cases are mainly 
characterized by the formation of extended plateaux cen- 
tered around the negative values. This result means that 
although the cross-correlation coefficient (fT3|) is close to 
zero for strongly magnetized flows (see Figure [3]), a wide 
range a values is often reached locally. 



2. Shear- and pseudo-Alfven wave decomposition 

In the presence of an external magnetic field, it is con- 
venient to describe the flow dynamics in terms of shear- 
and pseudo-Alfven waves, or in other words to use, re- 
spectively, the toroidal and poloidal components of the 
z ± fields (see equation (|SD). Indeed, the Alfven waves 
dynamics for the stronger magnetized flows have crucial 
consequences on the turbulent properties. We will use 
here the shear- and pseudo-Alfven wave decomposition 
to analyze our numerical simulations and, therefore, we 
will not considered the Bq — case anymore. 

In Figure [6l we show the temporal evolutions of ener- 
gies E± and associated, respectively, to the shear- 
Alfven and pseudo-Alfven waves; they are defined as 



(19) 



and are not inviscid invariants. Note that E ± ^ E^+E^ 1 
because the energy contained into the k± = modes 
are not included in the toroidal/poloidal decomposition 
(although it is, of course, in the original cartesian fields). 
First, we observe a slowdown of the energy decay when 
the intensity of Bq increases. It is a behavior similar to 
the one found in Figure [1] for the energies E ± . With such 
a decomposition, a similar behavior is also found for runs 
Ilia and IVa, and runs Va and Via. The important 
new information is about the initial increase of energies 
which is more pronounced for runs Va and Via, and 
for the shear- Alfven waves. These energies are in fact 
pumped from the k± = mode (the total energy is not 
an increasing function). Note that the same behavior is 
found for the — polarity. 

Figures [3 presents the temporal evolution of the dissi- 
pations 



^+ 2 (t) = z.<[Vx(z+ 2 )] 2 >(t), 



(20) 



for, respectively, the shear and pseudo-Alfven waves 
(only the + polarity is shown since the same behavior 
is found for the — polarity). No clear difference is found 



between the type of dissipation. We also note no sig- 
nificant difference with Figure [2] except a factor two in 
magnitude because here we do not see the total dissi- 
pation for a given polarity but either the shear- or the 
pseudo-Alfven waves contribution. 

Figure [8] presents the temporal evolution of the 
alfvenicity (or Alfven ratio) 



rA lt2 (i) 



(21) 



for the shear- and pseudo-Afven waves respectively, with 



Eh 



= 2 < «2 



> 



and 



E h2 



2 < ( < 2 



& 1,2 



) 2 > 



(22) 



(23) 



This plot is particularly interesting since it shows that 
shear- Alfven waves and pseudo-Alfven waves behave dif- 
ferently with an Alfven ratio of about one for the latter 
and significantly smaller than one for the former. Since 
for strongly magnetized flows the perpendicular fluctu- 
ations are mainly made of shear- Alfven waves and the 
parallel one made of pseudo-Alfven waves we have here 
a prediction that can be compared with measurements 
made in natural plasmas like in the solar wind. Addi- 
tionally, we observe the same oscillations as in Figure [H 
where the same type of analysis on the time-scales may 
be made. 

Figure [9] displays the temporal evolution of the spectral 
alfvenicity for shear- Alfven waves 



^A 1 ik h t) 



E\(k h ty 



(24) 



with fey = 0,1,2 and 3 (run Ila to IVa). The initial 
Afven ratio is close to unity for every parallel wavenum- 
bers, then a different behavior is found for the 2D state 
(kn = 0) which deviates strongly from the equipartition 
and tends approximately to 1 /2 independently of the Bo 
intensity. For the 3D modes, the spectral Alfven ra- 
tio oscillates around unity meaning a tendency towards 
equipartition between the kinetic and magnetic energies. 
This tendency is stronger for stronger magnetized flows. 
Thus the 3D modes follow the dynamics expected in wave 
turbulence in which an exact equipartition happens [3] . 
It is actually the 2D state which explains the behavior 
found previously in Figure [8] where a discrepancy from 
the equipartition was observed. Therefore, Figure [8] is 
not in contradiction with the wave turbulence regime and 
offer a new possible interpretation of observations in nat- 
ural plasmas like the solar wind. Note that the same type 
of results is found for pseudo-Alfven waves (not shown) 
with a deviation from the equipartition for the 2D state. 
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FIG. 5: Probability distribution functions of the cross-correlation (runs la to Via) initially (top left), at the maximum of the 
global dissipation (top right), when 77% of the total energy is dissipated (bottom left) and at the final time (bottom right). 



B. Characteristic length- and time-scales 

Figures \W\ and [TT] present the time evolutions of the 
perpendicular integral length-scales, defined as 



L+ = 

-1-1,2 



J Ei 2 (k±,k\\)/k±dk±dk\\ 

K2 

and the parallel integral length-scales 

/ E^ 2 (k±, fc||)/fc|| dk±dk\\ 



L+ 



E l,2 



(25) 



(26) 



for, respectively, the shear- and pseudo-Afven waves. We 
first note, for shear- Alfven waves, a decrease of the per- 



pendicular scales and an increase of the parallel one af- 
terward we observe a saturation. These behaviors may 
be interpreted as a direct cascade in the perpendicular di- 
rection and a possible inverse cascade in the parallel one. 
The saturation phase with length-scales approximately 
frozen means that the spectra are well developed. The 
case Bo = 1 deviates from this analysis because the mean 
field is not strong enough to impose a full anisotropic dy- 
namics; it can be compared with a previous study made 
for pure isotropic turbulence [34]. For pseudo- Alfven 
waves, the situation is less clear even if we still observe 
globally the same behavior as before in the initial phase. 
It is the saturation phase which is the most different with 
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FIG. 6: Temporal evolution of energies E± (top) and .E^" 
(bottom) of the shear- and pseudo-Alfven waves for runs Ila 
to Via. 



FIG. 7: Temporal evolution of the dissipations vQ^ (top) and 
isQ 2 (bottom) of shear- and pseudo-Alfven waves for runs Ila 
to Via. 



an apparent oscillation that can be related to the period 
found from the previous analysis made for Figure 0J 

Figure [12] presents the temporal evolution of the non- 
linear time 



NLi 



and the Alfven time 



-1-1,2 



L i,2 

B 



(27) 



(28) 



based on the shear- and pseudo-Afven waves dynamics. 
We mainly observe a decrease of the Alfven time when 
the strength of the uniform field Bo increases whereas 
the nonlinear time is not strongly affected. Note that 



the profiles of the nonlinear times for the case Bo = 5 
and Bo = 10 look similar with oscillations whose periods 
are approximately the same as before. This is simply due 
to the definition used to build the nonlinear time which 
includes the previous length-scales. It is also this defi- 
nition which explains the initial decrease (t < 2) of the 
nonlinear time since the perpendicular integral length- 
scales follow the same behavior. 

Figure [13] shows the temporal evolution (top-left) of 
the time-scales ratio 



xt(t) = 



NL 



(29) 



between the Alfven ([28]) and eddy turnover ([27]) times. 
Only the case of shear-Alfven waves is shown since the 
same behavior is found for pseudo-Alfven waves. This 
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FIG. 9: Temporal evolution of the spectral Alfven ratio for 
Q 5 ] 15 shear-Alfven waves at a given parallel wavenumber (runs Ila 

~ toIVa). 

time 



FIG. 8: Temporal evolution of Alfven ratio for shear- and 
pseudo- Alfven waves (runs Ila to IVa). 



new plots give a quantitative estimate of the balance be- 
tween the time-scales that we discussed in the introduc- 
tion. We clearly see that the balance is sub-critical (xt (t) 
stays well below unity) as the strength of Bo increases 
with a value which remains about constant during the 
time of the simulation. Figure [13] displays also the spec- 
tral ratio between the Alfven and nonlinear time-scales 
for shear-Alfven waves. It is defined as 



xi(k±,k\\)(t) 



k ± z (i (t) 



|| -DO 



(30) 



Zl (t) = y/Ei(k ± ,k\\)k ± k\\ . (31) 



with 



The previous definition (|29|) is based on a global estimate 
of the time-scales. This new definition is more precise 
since it allows to take into account the scale at which the 
times are denned. Then, each time evolution is associ- 
ated with a couple of (spectral) scales (k±,ku). Different 
couples have been tried and only those for which the ra- 
tio fcy) displays an extended plateau have been 
reported. It is basically for times between t = 2 and 
t = 4, a range of time during which the small-scales have 
been produced and the nonlinear interactions are still im- 



portant. Note that we still observe oscillations that can 
be explained in terms of Alfven times scales. 

In Figure [EJ we report each couple (k±,hu) and show 

the anisotropic scaling law fey ~ k 2 /_ 3 as a reference. 
We see that such law fits well the points which means 
that the sub-critical balance (observe again here with 
Xi(k±i &||) < 1) is ^ill well described by the anisotropic 
scaling law (pQ). This property may be understood by 
a heuristic model [26[ where the time-scale ratio x I s 
supposed to be constant at all scales but not necessarily 
equal to unity which allows to use the IK phenomenology 
instead of the Kolmogorov one [10]. (Note that the same 
behavior is found when x^(fcj_, fen) is considered.) 

The question of the validity of the anisotropic scaling 
law fcy ~ k 2 / 3 /Bo (we use here the formulation given in 
pr| which includes the uniform magnetic field) beyond 
the inertial range, and in particular at larger scales, may 
be addressed from these numerical simulations. A first 
answer is given in Figure [TH with the couple (k± = 4, fey = 
1) which is at the largest scales of the system but does 
not follow the anisotropic law. 



C. Generalized anisotropy angles 



To quantify the degree of anisotropy associated with 
the flow, we use the generalized Shebalin angles (see 0, 



11 




2 4 6 8 10 

time 

FIG. 10: Temporal evolution, up to t = 10, of perpendicular 
(top) and parallel (bottom) integral length- scales for shear- 
Alfven (+) waves; same runs as in Figure [7] 
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FIG. 11: Temporal evolution, up to t = 10, of perpendicular 
(top) and parallel (bottom) integral length- scales for pseudo- 
Alfven (+) waves; same runs and legend as in Figure [TOl 



and reference therein), denned as 

Efcj|q(M)l 2 
£*if|q(M)l 2 



tan a = 



(32) 



where q is a vector field, like v, b or z ± in Figure [15j 
We start initially with a 3D isotropic flow for which # q ~ 
54,74°. Figure [T5l shows that the temporal evolution of 
the different angles (for the different fields) is the same 
with a behavior depending mainly on the intensity of 
Bq. For Bo = the energy transfer is similar in all 
directions and the temporal evolution of Shebalin angles 
remains almost constant, close to its initial value. For 
Bo = 5 and Bo = 15, the Shebalin angles quickly increase 
and stabilize around 78°. Thus, and as expected, the 
anisotropy develops with Bo- However, the flow is not 
totally confined in planes perpendicular to Bo like for a 
purely bi-dimensional fluid for which the Shebalin angle 
is 90°. Note that a stronger anisotropy is produced for 
runs Va and Via (for which the Reynolds number is 
higher) with angles up to 83°. This is explained by the 
wider range of k± available for such runs. 

In Figure [161 we report the generalized Shebalin angles 
for the vector fields j, w and w ± , where w ± = V x z ± . 
The same behavior as before is found with apparently a 



slightly stronger anisotropy (with angles closer to 90°) for 
the highest values of Bo. This is explained by the fields 
used which are built on the rotational of the previous 
fields shown in Figure [15] and thus to a higher depen- 
dence of relation ([32]) in perpendicular wavenumbers (in 
k\ instead of k\). 



IV. SPECTRAL ANALYSIS 

A. Reduced spectra 

Figure [TT] (left) displays the one dimension (reduced) 
spectra 

E + (k x ) = J E+(k)dk y dk z , (33) 
E + (k y ) = J E+(k)dk x dk z , (34) 
E+(k z ) = J E + {k)dk x dk y , (35) 



for different Bo intensity, with the same initial condi- 
tion, and at times where the spectra are the most ex- 
tended (i.e. t ~ 2 for runs la and Ila, t ~ 3 for Ila 
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and IVa, t ~ 4 for Va and Via). It basically illustrates 
the different spectral transfers in the perpendicular and 
parallel directions when the strength of the uniform mag- 
netic field increases, whereas the x and y dependence is 
roughly the same. The equivalent spectra with polarity 
— is not shown since it gives the same picture. Note that 
the scaling at large-scales is not in contradiction with the 
initial condition discussed in Section Hi CI which concerns 
the modal spectrum. 



B. Energy fluxes 

In Figure [T71 (right) the associate reduced energy fluxes 
are given in k x , k y and k z . They are built from the 
cartesian z + fields. A constant flux is only found at the 
largest scales of the system. The presence of a negative 
flux is sometimes observed for a uniform field Bo > 5. 
This property may be linked to the increase of the parallel 
length-scale seen in Figure [TUJ In this case, the flux is 
clearly not constant which means that it is likely the 
result of a non-local interaction rather than an inverse 
cascade. Note that the same behavior is also found four 
the — polarity. 

The locality or nonlocality of the energy flux and trans- 
fer of runs la to IVa has been investigated recently [38[ 
by means of different geometrical wave number shells. 
It is shown that the interactions between the two coun- 
terpropagating Elsasser waves may become nonlocal for 
strong magnetized flows. In particular, the energy flux in 
the k± direction is mainly due to modes which interact 



with the plane fey = (with local interactions), while the 
weaker cascade in the parallel direction is due to modes 
which interact with ku = 1 (with possible nonlocal inter- 
actions) [H, [39j] . This property has been interpreted as 
a signature of a transition towards the weak turbulence 
regime during which the number of effective modes in the 
energy cascade is reduced. 

C. Anisotropic spectra 

Figure [18] shows anisotropic spectra for shear- Alfven 
waves (polarity +) at times for which turbulence is fully 
developed (t ~ 4). First, we see spectra Ef(k±) (top) 
which are defined as 

E?(k ± ) = J Ef(k^k\\)dk\\. (36) 

Then two other sets of spectra are given: Ef 1 (k±, fey = 0) 
and E^(k±,k\\ = 1) (middle and bottom panels respec- 
tively). The most interesting case seems to be the mid- 
dle panel, i.e. the spectra of the two-dimensional (2D) 
state, from which we see a clear inertial range where a 
power law may be extracted. An attempt is made to find 
this power law by computing the compensated spectra 
E+E~k m . Different values are proposed in the insets. 
We see that the 2D state is characterized by approxi- 
mately m = 14/3 which means on average a spectrum 
steeper than the Kolmogorov one. This scaling is clearly 
different from the value found for E^(k±) where the Kol- 
mogorov value m = 10/3 is better fitted. The last case 
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FIG. 13: Temporal evolution for shear- Alfven waves of the ratio between the Alfven and nonlinear time-scales based on the 
anisotropic IK phenomenology (top- left; for runs Ila to Via) and the spectral ratio between the same times (runs Ila to Via). 
Note the use of a logarithmic coordinate in the top-left and bottom-rigth panels. 
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do not exhibit significant differences with, for example, a 
wider inert ial range. In fact, the latter effect is easier seen 
for spectra plotted at fixed, but large, fey (fey > 1). Fi- 
nally note the difference between these spectra and those 
found in Figure [TT] (left) with an inertial range easier to 
determine in Figure [18] which may be attributed to the 
choice of the representation (anisotropic spectra instead 
of reduced spectra). 
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00 



D. Anisotropic scaling laws 



FIG. 14: Couples of points extracted from Figure [13] well 
fitted by the anisotropic scaling law k\\ ~ h±. 

E^(k± : ku = 1) is the most difficult one to analyze and 
no clear scaling appears. Note that the hyperviscous runs 



In order to extract a scaling law between parallel and 
perpendicular wavenumbers, we plot the modes 
corresponding to the equality E^{k±_) = E^(k\\) with 

Ef(k\\)= [ Ef(k ±: k\\)dk ± . (37) 
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FIG. 15: Temporal evolution of generalized Shebalin angles 
for the velocity, magnetic, and z fields (runs la to Via). 



FIG. 16: Temporal evolution of generalized Shebalin angles 
for the vorticity, current density and uj vorticity fields (runs 
la to Via). 



The result is given in Figure [19] for runs Ila (t ~ 2), 
Illa-IVa (t ~ 3) and Via (t ~ 4). We clearly see 
different slopes for different magnitudes of £?o, with an 
isotropic law fey ~ kj_ for Bo = 1 and an anisotropic law 

around fen ~ /c ± for £?o > 15 (see insets). For all cases, 
we see that the scaling law extends to the dissipative 
range. The same behavior is found for the pseudo-Alfven 
waves (not shown here) . Note that with this method the 
scaling law extracted suffers from an average effect since 
each spectrum is obtained after summation over the par- 
allel or the perpendicular direction. Nevertheless, the 
anisotropic prediction proposed by [10] is often recovered, 
but as it was explained above for a sub-critical balance 
between the Alfven and nonlinear times, which may be 
understood in a wider context [26[ as discussed in the 
introduction. 



V. VISUALIZATIONS 

A. Spectral space 

At a time at which energy spectra are fully developped, 
Figure [20] displays perpendicular, at fcy = 0, and paral- 
lel, at k y = 0, cuts in Fourier space for £? + (k) in flows 
at v = r] = 4.10" 3 without (B = 0, t = 2; run la) or 
with (Bo = 15, t = 3; run IVa) an applied magnetic 
field. Initially, for both flows, the isotropic energy injec- 
tion corresponds to spherical shells with maximum radius 



k = 8. The spectra then evolves depending on the level 
of the flow magnetization. Indeed, at Bo = 0, the max- 
imum spectral radius increases in all directions, mean- 
ing an isotropic energy transfer towards small scales, 
while at Bo = 15, the three-dimensional energy spectrum 
collapses into ellipsoidal shapes with ratio l/6th, corre- 
sponding to an anisotropic transfer, strongly inhibited in 
the Bq parallel direction. In this case, in Bq perpendic- 
ular planes (shown here at fen = 0), one can observe a 
loss of excitation at higher modes together with a loss of 
axisymmetry, with two preferred directions, compared to 
non-magnetized flows. 

Figure [2T1 shows the case of strongly magnetized flows, 
Bo = 30, at lower viscosity v = 10 -3 , and resolved with 
512 2 x 64 grid points (t = 4; run Via). The aspect ratio 
of the spectral ellipsoidal shape decreases up to l/10th 
and in transverse planes, a star shape with several "jets" 
appears. As time evolves (not shown), the number of 
these jets increases leading to an enhanced isotropy in 
tranverse planes (at kn > 0). In all flows, similar obser- 
vations stand for E~ (k) spectra. 



B. Physical space 

In order to understand the observed spectral structures 
in transverse planes for magnetized flows, we first visual- 
ize their spatial counterparts, once some Fourier ampli- 
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FIG. 17: Reduced spectra (left) E + built from the cartesian fields z + and the associated reduced energy fluxes (right) II + for 
the variables k x (solid line), k y (dashed line) and k z (dotted line for E + and long-dashed line for II + ). In the latter case when 
a negative flux is found, the absolute value is taken (dotted line). (Runs la to IVa and run Via, from top to bottom.) 



tudes at wavevectors (k Xl k yi k z = fen) are filtered for a 
given field. Hence structures only corresponding to the 
2D state are obtained with (k Xl k yi k\\ > 0) modes filtered, 
and structures for 3D modes (k\\ > 0) are obtained with 
(k x , k y , k\\ = 0) modes filtered. Figure [22l for a flow with 
Bq = 30 (run Via), displays vorticity and current iso- 



surfaces for the 2D state at the same time as Figure [2TJ 
t = 4, and at a later time t = 6. The transverse spec- 
tral star shape is related to the spatial distribution of the 
vorticity and current sheets, perpendicularly to two pecu- 
liar directions at t = 4, and more irregularly distributed 
at t = 6 for which higher number of jets is observed in 
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FIG. 18: Energy spectra of shear- Alfven waves E~± (solid line) and (dashed line) vs. the perpendicular wavenumbers k± 
after integration over all parallel wavenumbers (top), for k\\ = (middle) and for k\\ = 1 (bottom). Each of these three columns 
presents, from left to right respectively, runs Va and Via (y = 10 -3 with Bo = 15 and 30), and the hyperviscous run Vila 
(y = 10 -6 with Bq = 15). Inset: compensated product of energy spectra, E + E~k rn for a given m. 



spectral transverse planes (not shown). 

Similarly, for a flow with Bq = 15 (run IVa), vorticity 
and current isosurfaces for states with fen = and fen > 
are shown in Figure [23] at t = 7, when the total energy 
loss is about 40%. The 2D state structures are again re- 
lated to the star shape in transverse spectral planes (see 
Figure [20]h while the vorticity and current sheets with 
fen > present filamentary structures. From our knowl- 
edge it is the first time that such filaments are reported 
for (strongly) magnetized flows. 

When looking at the dynamics in physical space (with- 
out filtering) , shown in Figure [23J the vorticity and cur- 
rent intensities are superimposed sheet-like structures 
aligned along the ambiant magnetic field. At t = 7, a 
filament formation is observed within the sheets. This 
can be related to the filamentary structures with ku > 
(see Figure [23]) that do not exist in the 2D state, meaning 
that this sheet filamentation is mainly due to the wave 



components. At later times, t = 10, with a total energy 
loss of about 55%, the vorticity and current sheets are 
disruped by dissipation effects. 



Figure [25] shows the distribution of the cross- 
correlation values at z = tt (top) and in the entire nu- 
merical box (bottom). A comparison with the current 
and vorticity distribution shows that the high (absolute) 
value of the cross-correlation coincide with the position of 
dissipative structures which means that the velocity and 
magnetic field fluctuations tend to be aligned at small 
scales. This result corroborates recent works on the dy- 
namic alignment in MHD [40] where a statistical model 
is proposed. 



17 



100 



X 10 



10 



run256 B„=l 




100 



X 10 



run256 B = 5 



100 




10 



100 



100 



a 10 



run256 B =15 




10 



iiiiiinnffl - 
lll III 1 1 1 1 1 1 1 1 1 1 1 1 1 1 
llllllll it III III 1 1 1 1 

1 1 1 1 1 1 1 ir 

ill II Ml MM 



100 



X 10 



100 



run512 2 x64 B =30 



ki 



o ao ioo ibo aoo 



10 



iiiiiiiiimP "'" 

IIIIIIIIIHI 



100 
ki 



000 



FIG. 19: Anisotropic scaling laws between wavenumbers k± and k\\ (see text) for runs Ila to IVa and run Via. The inset 
displays two compensated scalings: k\\ (k±) k^ 1 (diamonds) and k\\ (k±) k^ 2 ^ 3 (crosses). 



VI. RUNS lb AND lib 

In a last set of simulations we change the initial condi- 
tion, as explained in Section Hi C 3[ to evaluate in partic- 
ular their influence on the dynamics. It is thought that 
this new initial condition is more appropriate to turbu- 
lent flows with a modal spectrum at large-scales (larger 
than the integral length-scale) in k\ which is in agree- 
ment with the phenomenology for freely decaying turbu- 
lence [22] . These runs correspond to a strong magnetized 
flow (Bo = 15) and high resolution (512 2 x 64). We will 
not focus on the temporal decay which has been ana- 
lyzed recently [22[ and we will only look at the spectral 
behavior. 

In Figure [26] we show the ID spectra E^ 2 (k±) f° r 
shear- and pseudo-Alfven waves integrated over all paral- 
lel wavenumbers. The time chosen is the one for which we 
have a fully developed turbulence. Despite the high res- 
olution no clear inert ial range appears. The Kolmogorov 
scaling is given as a reference that is roughly followed. 

The next Figure E3 gives at the same time the energy 
spectra E^ 2 and E± 2 of shear- and pseudo-Alfven waves 



for the 2D state (fcy = 0). The most remarkable result is 
the presence of a relatively extended inertial range char- 
acterized by a compensated energy spectrum on average 
around the IK prediction, i.e. in kj 3 ^ 2 (note however 
the presence of a bottleneck effect for the viscous run 
(left)). We conclude that the integration over the par- 
allel wavenumbers tends to hide the true scaling by an 
average effect. 



The last Figure [28] gives at the same time the en- 
ergy spectra E± 2 and E± 2 of shear- and pseudo-Alfven 
waves at a fixed parallel wavenumber (fey = 1). Once 
again a relatively extended inertial range is found. It is 
characterized by a compensated energy spectrum steeper 
than the previous one with an index around kj 2 and 

kj 7 ^ 3 for respectively the hyperviscous and viscous case. 
Other spectra at higher fixed parallel wavenumbers are 
not shown because there are characterized by a smaller 
inertial range from which it is difficult to find a power 
law scaling. 
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FIG. 20: E + (k) cuts in Fourier space at k\\ — (left) and k y — (right) for flows at Bo — (top; run la) at t = 2, and Bo = 15 
(bottom; run IVa) at t = 3. Color bars normalized to 1 for the maximum intensity (red) and to for the minimum (blue) one. 




FIG. 21: E + (k) cuts in Fourier space at k\\ = (left) and k y = (right) for flows at B = 30 (run Via using 512 2 x 64 grid 
points) at t — 4. Color bars normalized to 1 for the maximum intensity (red) and to for the minimum (blue) one. 

VII. SUMMARY AND CONCLUSION magnetic field B is investigated. A parametric study in 

terms of Bo intensity is made to show the development of 
anisotropy. In general, the temporal evolutions show os- 
In this paper, we present a set of 3D direct numer- dilations that are associated with the presence of Alfven 
ical simulations of incompressible decaying MHD tur- waves. The dynamics is slower for strongly magnetized 
bulence in which the influence of an external uniform 
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B =30, t=4 



B =30, t=6 




FIG. 22: Isosurfaces of vorticity (blue) and current (orange) intensities for the 2D state k\\ = (see text) for a flow with 
Bo = 30 (run Via), drawn at 20% of their respective maxima. At t = 4 (left), |w| macc = 8.8 and \j\max — 11.3, and at t = 6 
(right) \w\max = 6.6 and |j| macc = 8.1. 



B =15, t=7 



B =15, t=7 




FIG. 23: Filtered vorticity (blue) and current (orange) intensities with k\\ — (2D state, left) and k\\ > (right) for Bo — 15 
(run IVa) : isosurfaces are drawn at 27% and 20% of their respective maxima (/cm = 0; |w| macc = 5.7 and |j| m acc = 9., and 
k\\ > 0; \-w\max = 18.5 and |j| maaJ = 19.2). 



flows with, in particular, a cross-correlation between the 
velocity and the magnetic field fluctuations frozen on av- 
erage around its initial (small) value but with, locally, 
a wide range of possible values. For all temporal re- 
sults, one can see that the flows with the highest values 



of (> 5) behave quite similarly while for Bq = 1, 
the flow presents a transient regime between the case 
without background magnetic field and the other cases. 
We also discuss the presence of a sub-critical balance be- 
tween the Alfven and nonlinear times with both a global 
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FIG. 24: Temporal evolution of vorticity (blue) and current (orange) intensities for the same run as in Figure [23] Isosurfaces 
drawn at 27% and 20% of |w| macc and \j\max instantaneous maxima. At t = 0; |w| macc = 2.7 and \j\ ma x = 3., at t = 4; 
|w| macc - 24.3 and \j\max = 27.6, at t = 7, |w| macc - 17.9 and \j\max = 25.6, and at t = 10, |w| macc - 9.2 and |j| moa : - 15.4. 



and a spectral definition. This regime is still associated 
with the anisotropic scaling laws (pQ) between the per- 
pendicular and the parallel wavenumbers. The nonlinear 
dynamics of strongly magnetized flows is characterized 
by a different k± -spectrum if it is plotted at a fixed fey 
(2D spectrum) or if it is integrated (averaged) over all fen 
(ID spectrum). In the former case a much wider iner- 
tial range is found with a steep power law, closer to the 
wave turbulence prediction than the Kolmogorov one like 
in the latter case. Note that the inertial range of these 



spectra is better seen for the shear- and pseudo-Alfven 
waves rather than for the cartesian fields. 

One of the most important results of this paper is the 
difference found between the fcx-spectra plotted after in- 
tegration over fen and those at a given fen. This point 
is generally not discussed in numerical works whereas it 
appears to be a fundamental aspect of this problem. Di- 
rect numerical simulations of the Alfven wave turbulence 
regime seems to be still out of the current numerical ca- 
pacity [32[ and only the detection of the transition to- 
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FIG. 25: Top-left: isosurfaces (at z — tt) of the cross-correlation coefficient; Top-right: the corresponding isosurfaces of the 
current (red) and vorticity (blue) drawn at 21% of their respective maxima; Bottom- left: isosurfaces of the cross-correlation 
coefficient at —0.70 (yellow), —0.75 (blue) and —0.90 (red); Bottom-right: isosurfaces of the cross-correlation coefficient at 0.96 
(yellow) and -0.96 (pink). (Run IVa at time t = 10.) 



wards such a regime seems possible. In such a study it 
is crucial to avoid any noisy effect linked, for example, 
to the initial condition (or forcing) that could favored 
one particular type of spectrum. But the other effect 
that could hide the true dynamics of strongly magnetized 
flows is the averaging effect as we have clearly seen in 
the second set of simulations: the presence of an inertial 
range was not obvious from a first global analysis (Fig- 
ure [26]) whereas it was clear from 2D spectra (Figure [27]). 



This averaging effect may be due to the moderate spatial 
resolution used but also to the regime which is in a tran- 
sition phase towards the wave turbulence regime. If we 
extrapolate such a result to natural plasmas like the one 
found in the interplanetary medium (inner solar wind) 
then we may interpret the current spectra as averaging 
spectra (since we are not able to report spectra at a given 
parallel wavenumber with only one spacecraft). Then it 
is not surprising that we observe both an anisotropic flow 
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k ± 

FIG. 26: Energy spectra E^ 2 (k±) for shear- (solid) and 
pseudo- Alfven (dashed line) waves integrated over all paral- 
lel wavenumbers in the viscous case (lb). The straight line 
follows a k^ 5 ^ 3 law. 

with approximately a Kolmogorov scaling. 

In a recent numerical analysis [33[ dedicated to the de- 
velopment of anisotropy and wave turbulence in forced 
incompressible reduced MHD flows, a change of spectral 
slope was reported for the /c^-energy spectrum when the 
forcing is applied on a larger range of parallel wavenum- 
bers with no driving of the fen =0 modes. In the light 
of the present paper this finding may be interpreted as 
a way to decrease the averaging effect which is mainly 
due to the dissipative scales. Indeed, when a larger 
parallel wavenumbers is excited the spectrum integrated 
over all fey is more sensitive to the non dissipative par- 
allel wavenumbers and tends therefore to reveal the true 
scaling. Another way to avoid this noisy effect would 
have been to plot the spectra at a given but low parallel 
wavenumber in order to avoid the dissipative range. 

The question of the power law index predicted by the 
wave turbulence theory has not been addressed so far. A 
/cj 2 -spectrum is expected for strongly magnetized flows 
in the regime of wave turbulence (even without assuming 
a restriction on k± and fey |41|). In our simulations a 
scaling close to this value is found when hyperviscosity 
is used (Figure [28]) in the second set of simulations at 
fcy = 1 whereas the 2D state (spectrum at fcy = 0) scales 

3/2 

on average around k ± 1 . This latter result is the same 
as the one found generally in 2D isotropic MHD turbu- 
lence HHH. The steep power law reported in k ^ with 
a £ [2, 7/3] may be attributed to the very first sign of the 



wave turbulence regime that should be confirmed never- 
theless at higher resolution. The a = 7/3 case is a priori 
unexpected although it was seen as a transient re gim e 
before the finite energy flux solution settles down [14{ . 
However, no change of slope is observed in our simula- 
tion because, in particular, for larger times the reduction 
of the inertial range does not allow to conclude about the 
inertial scaling law. The case a = 7/3 is also a solution 
predicted by a heuristic model based on a sub-critical 
balance between the Alfven and the nonlinear times [26[ . 
In this case, the total energy spectrum satisfies the rela- 
tions E(k±, fey) ~ kj_ a k7 p , with 3a + 2/3 = 7. Thus the 
a = 7/3 solution implies no fen -scale dependence which 
could be linked to the weakness of parallel transfers. 

Our analysis in the physical space has revealed an im- 
portant new information about structures. A filament 
formation is observed within the current and vorticity 
sheets. This important property may be explained by 
the specific condition of our simulation (large Bo and 
large Reynolds number) and has to be confirmed at 
higher Reynolds numbers. The classical picture of cur- 
rent sheets in MHD turbulence may be wrong in the 
strongly anisotropic case and filaments are may be the 
right picture. This result may be compared with astro- 
physical plasmas like in the solar corona where extremely 
thin (dissipative) coronal loops (filaments or "strands") 
are observed. Although their presence is well accepted, 
the origin of these filaments is still not well explained. 
Turbulence and Alfven wave could be the main ingredi- 
ents m« 

Other questions about scaling laws for structure func- 
tions and intermittency for strongly magnetized flows are 
not discussed here. Forcing numerical simulations are 
then necessary which is out of the scoop of this paper. 
The unbalance case has not been addressed in this paper. 
It is also an important issue not only from a theoretical 
point of view but also from an observational point of 
view since the most analyzed astrophysical plasma, the 
inner solar wind, is mainly made of outwards propagating 
Alfven waves. This point is left for future works. 



Acknowledgments 

We thank A. Alexakis for useful discussions. This 
work is supported by INSU/PNST-PCMI Programs and 
CNRS/GdR Dynamo. This work was supported by 
the ANR project no. 06-BLAN-0363-01 "HiSpeedPIV" . 
Computations time was provided by IDRIS (CNRS) 
Grant No. 070597. 



[1] T. Tajima, and K. Shibata, Plasma Astrophysics (West- (1999). 

view Press, Boulder, USA 2002). [3] S. Galtier, J. Low Temp. Phys. 145, 59 (2006). 

[2] M.L. Goldstein and D.A. Roberts, Phys. Plasmas 6, 4154 [4] PS. Iroshnikov, Soviet Astron. 7, 566 (1964). 



23 



run512 2 X64 B =15 (lb) run512 2 X64 B =15 (lib) 




1000 



1 10 100 1000 1 10 100 1000 

FIG. 27: Energy spectra E^ 2 (solid) and E± 2 (dashed) of shear- (top) and pseudo- (bottom) Alfven waves for the 2D state 
(k\\ — 0). The viscous (lb) (left) and the hyperviscous case (lib) (right) are shown. Inset: compensated energy spectra 
El 2 (k ± ,0)E- 2 (k ± ,0)k m . 



[5] R.H. Kraichnan, Phys. Fluids 8, 1385 (1965). 
[6] D. Montgomery, and L. Turner, Phys. Fluids 24, 825 
(1981). 

[7] J.V. Shebalin, W.H. Matthaeus and D. Montgomery, J. 

Plasma Phys. 29, 525 (1983). 
[8] J.-C. Higdon, Astrophys. J. 285, 109 (1984). 
[9] S. Oughton, E.R. Priest and W.H. Mattaheus, J. Fluid 

Mech. 280, 95 (1994). 
[10] P. Goldreich, and S. Sridhar, Astrophys. J. 438, 763 

(1995) . 

[11] C.S. Ng and A. Bhattacharjee, Astrophys. J. 465, 845 

(1996) . 

[12] R.M. Kinney and J.C. McWilliams, Phys. Rev. E 57, 
7111 (1998). 

[13] W.H. Matthaeus, S. Oughton, S. Ghosh and M. Hossain, 

Phys. Rev. Lett. 81, 2056 (1998). 
[14] S. Galtier, S.V. Nazarenko, A.C. Newell and A. Pouquet, 

J. Plasma Phys. 63, 447 (2000). 
[15] S. Galtier, S.V. Nazarenko, A.C Newell and A. Pouquet, 

Astrophys. J. 564, L49 (2002). 



[16] L.J. Milano, W.H. Matthaeus, P. Dmitruk and D.C. 
Montgomery, Phys. Plasmas, 8, 2673 (2001). 

[17] W.-C. Muller, D. Biskamp and R. Grappin, Phys. Rev. 
E 67, 066303 (2003). 

[18] M.K. Verma, Phys. Rep. 401, 229 (2004). 

[19] B.D.G. Chandran, Phys. Rev. Lett. 95, 265004 (2005). 

[20] W.-C. Muller and R. Grappin, Phys. Rev. Lett. 95, 
114502 (2005). 

[40] S. Boldyrev, Phys. Rev. Lett. 96, 115002 (2006). 

[22] B. Bigot, S. Galtier and H. Politano, Phys. Rev. Lett. 
100, 074502 (2008). 

[23] T.S. Horbury, in Plasma Turbulence and Energetic Par- 
ticles in Astrophysics, edited by M. Ostrowski and R. 
Schlickeiser (U. Jagiellonski, Cracow, 1999) p. 115. 

[24] S. Dasso, L.J. Milano, W.H. Matthaeus and C.W. Smith, 
Astrophys. J. 635, L181 (2005). 

[25] B.G. Elmegreen and J. Scalo, Annu. Rev. Astron. As- 
trophys. 42, 211 (2004); J. Scalo and B.G. Elmegreen, 
Annu. Rev. Astron. Astrophys. 42, 275 (2004). 

[26] S. Galtier, A. Pouquet and A. Mangeney, Phys. Plasmas 



24 



run512 2 X64 B =15 (lb) 



+ 




run512 2 X64 B =15 (lib) 



^4 



+ 



10" 



10" 



10" 



10" 



-10 



10 
IQ" 1 



10 100 
run512 2 X64 B =15 (lb) 




10 100 



1000 



1000 



b< 

b? 
b? 

+ 




b? 

+ 

b^ 
+ 

CV2 



10" 



10" 



10" 



10" 



10 



10 



-10 



-12 



10 100 
run512 2 X64 B =15 (lib) 




1000 



1000 



FIG. 28: Energy spectra i?^ 2 (solid) and E 12 (dashed) of shear- (top) and pseudo- (bottom) Alfven waves for k\\ = 1. The vis- 
cous (lb) (left) and the hyperviscous case (lib) (right) are shown. Inset: compensated energy spectra E^ 2 (k±, l)E^ 2 (k±, l)k m . 



12, 092310 (2005). 
[27] J. Cho and E. T. Vishniac, Astrophys. J. 539, 273 (2000). 
[28] J. Maron and P. Goldreich, Astrophys. J. 554, 1175 

(2001). 

[29] D. Shaikh and G. Zank, Astrophys. J. 656, L17 (2007). 
[30] V.E. Zakharov, V. L'vov and G.E. Falkovich, Kolmogorov 

Spectra of Turbulence I: Wave Turbulence. (Springer- 

Verlag, Berlin, Germany 1992). 
[31] A.C. Newell, S.V. Nazarenko and L. Biven, Physica D 

152-153, 520 (2001). 
[32] S.V. Nazarenko, New J. Phys. 9, 307 (2007). 
[33] J.-C. Perez and S. Boldyrev, Astrophys. J. 672, L61 

(2008). 

[34] S. Galtier, H. Politano and A. Pouquet, J. Plasma Phys. 
61, 507 (1999). 

[35] S. Oughton, P. Dmitruk, and W.H. Matthaeus, Phys. 
Plasmas 11, 2214 (2004). 



[36] L. Smith, and F. Waleffe, Phys. Fluids 11, 1608 (1999). 
[37] R. Bruno, and V. Carbone, Living Rev. Solar Phys. 2, 1 
(2005). 

[38] A. Alexakis, B. Bigot, H. Politano and S. Galtier, Phys. 

Rev. E 76, 056313 (2007). 
[39] A. Alexakis, Astrophys. J. 667, L93 (2007). 
[40] J. Mason, F. Cattaneo and S. Boldyrev, Phys. Rev. Lett. 

97, 255002 (2006). 
[41] S. Galtier and B.D.G. Chandran, Phys. Plasmas 13, 

114505 (2006). 

[42] H. Politano, A. Pouquet and V. Carbone, Europhys. Lett. 
43, 516 (1998). 

[43] D. Biskamp and E. Schwarz, Phys. Plasmas 8, 3282 
(2001). 

[44] B. Bigot, S. Galtier and H. Politano, accepted for Astron. 
Astrophys. 



